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Abstract: The signals expected in WIMP direct detection experiments depend on the 
ultra-local dark matter distribution. Observations probe the local density, circular speed 
and escape speed, while simulations find velocity distributions that deviate significantly 
from the standard Maxwellian distribution. We calculate the energy, time and direction 
dependence of the event rate for a range of velocity distributions motivated by recent 
observations and simulations, and also investigate the uncertainty in the determination 
of WIMP parameters. The dominant uncertainties are the systematic error in the local 
circular speed and whether or not the MW has a high density dark disc. In both cases 
there are substantial changes in the mean differential event rate and the annual modulation 
signal, and hence exclusion limits and determinations of the WIMP mass. The uncertainty 
in the shape of the halo velocity distribution is less important, however it leads to a ~ 5% 
systematic error in the WIMP mass. The detailed direction dependence of the event rate 
is sensitive to the velocity distribution. However the numbers of events required to detect 
anisotropy and confirm the median recoil direction do not change substantially. 
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1. Introduction 

Weakly Interacting Massive Particles (WIMPs) are a well motivated dark matter candi- 
date |Q, ||, ||]. They can be directly detected in the lab via their elastic scattering off 
target nuclei in dedicated detectors § and experiments are now probing the theoretically 
favoured regions of parameter space ||, 0] . 

The nuclear recoil event rate is energy, time and direction dependent. Due to the 
Earth's orbit about the Sun the net velocity of the lab with respect to the Galactic rest 
frame varies over the course of a year. The net speed is largest in the Summer and 
hence there are more high speed WIMPs, and less low speed WIMPs, in the lab frame. 
This produces an energy dependent annual modulation in the event rate with amplitude 
of order a few per-cent j8j. The WIMP flux in the lab frame is sharply peaked in the 
direction of motion of the Sun (towards the constellation CYGNUS), and hence the recoil 
spectrum is peaked in the direction opposite to this ||. This directional signal is far 
larger than the annual modulation; the event rate in the backward direction is roughly an 
order of magnitude larger than that in the forward direction |J. A detector which can 
measure the recoil directions is required to detect this signal (see Ref. [jnj] for an overview 
of the status of directional detection experiments). The time and direction dependence of 
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the event rate are signals which can be used to discriminate WIMP induced recoils from 
backgrounds, while the energy dependence of the event rate can be used to measure the 
WIMP mass gl], p], 0, 0, [D| |]| . 

H, time §Si| 



The energy 



2|, ||, ||, |§ and direction |27|, ||, |9] 



dependence of the event rate all depend on the ultra-local WIMP distribution. Historically 
analytic halo models have been used to calculate the WIMP signals and analyse data. In 
the past few years velocity distributions from high resolution simulations of the formation 
of Milky Way like halos |?C], 31, 32, 33, 34, [3f| have been used to calculate some of the 

m. 



WIMP signals pL 37, 38, B2, 39, 34, 35, 40 



In this paper we study the energy, time and direction dependence of the elastic scat- 
tering event rate for a range of velocity distributions motivated by recent observations and 
simulations. Where possible we use WIMP and target nuclei independent parameterisa- 
tions of the event rate. We also investigate the uncertainty in exclusion limits, determi- 
nations of the WIMP mass and the number of events required to detect the directional 
signal. For studies of the impact of uncertainties in the WIMP velocity distribution on 
the interpretation of data from recent experiments see Refs. (3^, |34|, |35], 41, 42] for 
elastic scattering, Refs. 34, [35|, ^] for inelastic scattering and Ref. [41] for momentum 
dependent scattering. 

In Sec. |2| we review the calculation of the event rate. In Sec. || we discuss recent results, 
in particular from numerical simulations, on the local dark matter distribution and present 
the velocity distributions which we will use to calculate the WIMP signals in Sec. ||[ Finally 
we conclude with a summary in Sec. |. 



2. Event rate 

Assuming spin-independent coupling the differential event rate (number of events per unit 
energy, per unit time, per unit detector mass) is given by [|l], |llj 

§ ( E,t ) = -^AV(E ) r (2.1) 
etc/ zfi px m x Jv min v 

where p x is the ultra-local WIMP density, f(v,t) the normalised ultra-local WIMP speed 
distribution in the rest frame of the detector, <t p the WIMP scattering cross section on the 
proton, /i px = (m p m x ) / (m p + m x ) the WIMP-proton reduced mass, A and F(E) the mass 
number and form factor of the target nuclei respectively and E is the recoil energy. The 
lower limit of the integral, v m { n , is the minimum WIMP speed that can cause a recoil of 
energy E: 

where txia is the atomic mass of the detector nuclei and \xp^ x the WIMP-nucleon reduced 
mass. 
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The WIMP speed distribution in the detector rest frame is calculated by carrying out, 
a time dependent, Galilean transformation: v — > v = v + v e (i) . The Earth's motion 
relative to the Galactic rest frame, v e (t), is made up of three components: the motion of 
the Local Standard of Rest (LSR), vlsr, the Sun's peculiar motion with respect to the 
LSR, Vq, and the Earth's orbit about the Sun, v° rb . We use vlsr = (0, u c ,0), where v c 
is the local circular speed, and Vq = (11.1, 12.2, 7.3) kms -1 4S] in Galactic co-ordinates 



(where x points towards the Galactic center, y is the direction of Galactic rotation and z 
towards the North Galactic Pole) and the expressions for the Earth's orbit from Ref. 

The differential event rate, eq. (|2.l|) , depends on the target nuclei mass and also the (a 
priori unknown) WIMP mass. It is therefore useful (c.f. Ref. |l7]]) to consider the 'WIMP 
and target independent' quantity 

T(v min ,t) = (220 km s" 1 ) f°° dv . (2.3) 



The prefactor here is chosen to make r(u m i n , t) dimensionless (and of order unity) while 
allowing variations in the local circular speed (see Sec. [3.1.2 ). The differential event rate 
can then be written as 

where the prefactor 

C( X ,A) = -^A 2 F 2 (E), (2.5) 

contains the WIMP and target dependent terms and is independent of the astrophysical 
WIMP distribution. 

The direction dependence 0] of the event rate is most compactly written in terms of 



the radon transform of the WIMP velocity distribution [47] 

a^o = ^i^F\E)f(v mm , q) , (2.6) 

where dfi = d0dcos7, q is the recoil direction and /(« m m,q) is the 3-dimensional Radon 
transform of the WIMP velocity distribution /(v) 

f(v min , q) = J <5(v.q ~ %i„)/(v)d 3 « . (2.7) 

Geometrically the Radon transform, /(u m i n ,q), is the integral of the function /(v) on a 
plane orthogonal to the direction q at a distance v m i n from the origin. See Ref. (2^] for an 
alternative, but equivalent, expression. 

While the directional recoil rate depends on both of the angles which specify a given 
direction, the strongest signal is the differential of the event rate with respect to the angle 
between the recoil and the direction of solar motion, 7, || 

d 2 R P X °? >■-* 



dE d cos 7 47T/z px m x 



A F Z {E) / f(v min ,e0d4>, (2.8) 



1 Formally the effects of gravitational focusing by the Sun should be taken into account |44), however 
the resulting modulation in the differential event rate is small and only detectable with a very large number 
of events ^] . 
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where <j> is the azimuthal angle. As with the non-directional differential event it is useful 
to consider a dimensionless 'WIMP and target independent' quantity 

,-2-K 

T(v min , cos 7) = (220 km s" 1 ) / /> min , q) d<^> , (2.9) 



(2.10) 



so that eq. ( |2.8| ) can be rewritten as 

d 2 R C( X ,A)p x 



d£dcos7 2vr(220kms" 
where C(x, A) is defined in eq. 



jr T(v m in,cos7) , 



3. Dark matter distribution 

The traditional benchmark model for direct detection event rate calculations is the standard 
halo model (an isotropic isothermal sphere with p(r) oc r~ 2 ) which has velocity distribution 



/(v) oc < 



where a is the velocity dispersion. 



exp 




M 2 
2^ 



if |v| < v esc 
otherwise . 



(3.1) 



In Sec. 3.1 we review the status of observational measurements of the local dark matter 



density, circular speed and escape speed, and their relevance to the event rate calculation. 



In Sec. 3.2 we review results on the shape of the dark matter velocity distribution from 



recent numerical simulations. We conclude in Sec. 3.3 with a summary of the benchmark 
models which we will use to calculate the direct detection signals in Sec. ||. 

3.1 Observations 
3.1.1 Local density 

The local density appears in the normalisation of the event rate. We use the standard value, 
p x ~ 0.3 GeV cm -3 . For a given model for the MW it is possible to determine p x to ~ 10% 
accuracy [56, 57], however analyses which use a range of models or model independent 
methods find substantially larger, of order unity, uncertainties pi 



5£ 



] . Furthermore 

observations determine the local density averaged over a spherical shell, and the DM density 
in the stellar disc of simulated halos is ~ 20% larger than the shell average [ |BT| . 

Since the normalisation is proportional to the product of p x and a p , the uncertainty 
in p x translates directly into an uncertainty in constraints on, or measurements of, a p . As 
this uncertainty is the same for all experiments we do not explicitly consider it. 

3.1.2 Local circular speed 

The local circular speed, the speed with which stars orbit the Galactic centre at the solar 
radius, is related to the velocity dispersion by the Jean's equation [49|. 



Id(paf) 
P 



dr 



+ 2 



r 



(3.2) 
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where a r is the radial velocity dispersion and /3 = 1 — (cr 2 , + a^)/2a^ is the anisotropy 
parameter. For the standard halo model the velocity dispersion is isotropic (a r = ag = 
so that (5 = 0) and independent of radius and p(r) oc r -2 so that a r = v c j\[2. 

The standard value of v c is v c = 220kms _1 |50fl . A recent analysis using Galactic 
masers found a significantly higher value, v c = (254 ± 16)kms~ 1 , [51] and this has been 



adopted is some subsequent direct detection work p2|. It has been argued, however, that 
this analysis used overly restrictive models. Bovy et al. found v c = (236 ± ll)kms~ , 
assuming a flat rotation curve [p3[| , while McMillan and Binney find values ranging from 
v c = (200 ± 20)kms _1 to v c = (279 ± 33)kms _1 depending on the model used for the 
rotation curve |5^| . 

Given the significant systematic uncertainties in determinations of v c , we retain v c = 
220kms _1 as the default value, and also consider v c = 200 and 280 km s -1 . 

3.1.3 Local escape speed 

Particles with speed, in the Galactic rest frame greater than the local escape speed, v esc = 
\J2\Q(Rq)\ where <£(r) is the potential, are not gravitationally bound. The standard halo 
model formally extends to infinity and therefore the speed distribution has to be truncated 
at v esc 'by hand' (see e.g. Ref. |8|). Historically the standard value for the escape speed 
was v esc = 650 km s" 1 . A more recent analysis, using high velocity stars from the RAVE 
survey, finds 498 km s -1 < v eS c < 608 km s -1 with a median likelihood of 544 km s -1 [55]. 



We use as a default the median RAVE value, v esc = 544 km s , and also consider the 
old standard value, v esc = 650 km s -1 . 

3.2 Simulations 

High resolution dark matter only simulations of the formation of Milky Way like dark 
matter halos, find speed distributions which deviate systematically from a multivariate 



Gaussian (the simplest anisotropic generalisation of the Maxwellian distribution) |30| , 38 



32, 34]. There are more low speed particles, and the peak in the distribution is lower. 



The deviation is smaller in the lab frame than in the Milky Way rest frame however 



Kuhlen et al. [34| study the velocity distributions of the particles in the Via Lactea 2 
(VL2) |3j| and GHALO |34| simulations. In each case they consider the particles centered 
in a 1 kpc shell centered on galactocentric radius r = 8.5 kpc and also 100 sample spheres 
of radii 1 or 1.5 kpc, each centered at a point with r = 8.5 kpc. The shell contains a large 
number (~ 10 6 ) particles, allowing the average velocity distribution at the solar radius to be 
measured with small statistical errors. The spheres contain a smaller number of particles 
(~ 10 4 ), and hence have larger statistical errors, but are sensitive to local variations in 
the velocity distribution on ~ kpc scales. They find that the radial, v r and tangential, 
Vt = \/ v e + v % velocity distributions are well fit, apart from at large v t , by modified 



gaussian distributions [38]: 



/Or) = ^ ex P 
/K) = -|exp 



at 



(3.3) 
(3.4) 
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where N T / t are normalisation factors. The velocity distributions also have stochastic fea- 
tures at high speeds. There are broad bumps which vary from halo to halo, but are inde- 
pendent of position within a given halo and are thought to reflect the formation history of 
the halo Kuhlen et al. @ also find narrow spikes in some locations, corresponding 

to tidal streams. 

We use the shell and sphere median, 16th and 84th percentile fit parameters for the 
VL2 |33] simulation, given in Table 1 of Ref. As discussed by Kuhlen et al., the 

velocity dispersion and circular speed of these dark matter only simulations is lower than 
expected in reality; baryonic contraction will deepen the potential well and increase the 
circular speed. For VL2 the most likely speed is vq = 184 km s -1 and v c /vq 0.85. To 
allow a comparison with the standard halo model (with standard parameters) we scale the 
fit parameters so that the peak of the speed distribution matches that of the standard halo 
model. We truncate the fitted velocity distributions 'by hand' at the median value from 
RAVE, 

^vesc — 544 km s 1 . For the time averaged differential event rate we also consider 
the tabulated data from the VL2 simulation from Ref. [62]. As discussed by Ref. |P|] in 
this case the scaling results in some of the high speed features in the distribution being 
pushed beyond the escape speed. These features are likely to be significant for experiments 
which are only sensitive to the high speed tail of the speed distribution. 

It should be cautioned that the scales resolved by simulations are many orders of mag- 
nitude larger than those probed by direct detection experiments. Vogelsberger and White 
have developed a new technique to study the ultra- local dark matter distribution |55| . They 
find that the ultra-local dark matter distribution consists of a huge number of streams and 
is essentially smooth. Schneider et al. [36j have reached similar conclusions by studying 



the evolution of the first, Earth mass, microhalos to form |67L 68, p9[. This suggests (see 



also Ref. |70j ) that the ultra- local dark matter density and velocity distribution should not 
be drastically different (i.e. composed of a small number of streams) to those on the scales 
resolved by simulations. The ultra-local velocity distribution may, however, contain some 
fine-grained substructure jTl], ffi| . 

The simulations discussed above contain dark matter only, while baryons dominate in 
the inner regions of the Milky Way. Simulating baryonic physics is extremely difficult, and 
producing galaxies whose detailed properties match those of real galaxies is an outstanding 
challenge. Some recent simulations have found that late merging sub-halos are preferen- 
tially dragged towards the disc, where they are destroyed leading to the formation of a 
co-rotating dark disc (DD) |3l], 35]. 



The properties (density and velocity distribution) of the DD are highly uncertain. We 
consider 3 benchmark models, which aim to broadly span the range of plausible properties. 
The first benchmark model follows Ref. [37] modelling the DD velocity distribution as a 
gaussian with isotropic dispersion, odd = 50kms~ 1 and lag t>i ag = 50kms -1 , matching 
(roughly) the kinematics of the Milky Way's stellar thick disc. We use the central value of 
the DD density from Ref. |}7]], /?dd = Ph> where pu is the local halo density. 

Ref. [^3| argues that to be consistent with the observed morphological and kinematic 
properties of the Milky Way's thick disc, the Milky Way's merger history must be quiescent 
compared with typical ACDM merger histories. Hence the DD density must be relatively 
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label 


properties 


SHSP 


KtanHavH halo mnHpl fSTT^ witli 7? — 990 km ^ 7? — 554 km k 

DUCtllUClil U lilUUtl 1 Ull / VV lull £j £-i\j ±\.L1L O , CSC ' ' ) * 1V111. O 


SHt>cs C H 


SH with ?j c = 220km s~ , v csc = 650kms _1 


SHu T, 


with 7i — 900 km 7i — ^iVm 


SRv c R 


SH with v c = 280 km s -1 , v esc = 554 km s -1 


SIMsh 


scaled modified Maxwellian fit to VLz shell data, eq. (p,3[) irom Kel. 34): 
u r = 1.2 x 202 km s -1 , v t = 1.2 x 129 km s" 1 , a r = 0.93, a t = 0.64. 


bllvlspmea 


scaled median modified Maxwellian fit to VL2 sphere data: 

71 — 1 9 y 900 km « — 1 il - 1 9 Y 1 ^Ui a -1 rv — Q4 rti — fl fifi 




cpii oH l r^"t" n T^fiT^pTiti £>G iii r~i H 1 Ti ori ^/ltl "Vinrti 1 1 1 Q n irt" t \/ 1 ,v cp\ n pyp ri q "t" q ■ 
DCdlcCL 17C1 CCilLllco illULLlllcU. IVlclA. Wcllldll 11U L*J V ±JZ> a C iXOtbOt. 

v r = 1.2 x 186kms~\ v t = 1.2 x 124 km s" 1 , a r = 0.88, a t = 0.64. 


SIMsp84 


scaled 84th percentiles modified Maxwellian fit to VL2 sphere data: 
v r = 1.2 x 213 km s -1 , v t = 1.2 x 149 km s" 1 , a r = 0.99, a t = 0.67. 


DDpHcrL 


SHSP plus a dark disk with i>i ag = 50kms _1 , pdd = PR, cdd = 50kms _1 


DDpLcrL 


SHSP plus a dark disk with vi ag = 50kms _1 , Pdd = 0.15/9H, ct dd = 50km s -1 


DDpLo-H 


SHSP plus a dark disk with vi ag = 50kms _1 , pdd = 0.15/9H, odd = 100 kms -1 



Table 1: Summary of benchmark models. 



small, pdd < 0.2/?h, at the lower end of the range of values considered in Ref. |37]]. Refs. [73 



4C] also argue that the DD velocity dispersion is likely to be substantially larger than 



that of the stellar thick disc. In order to study the effects of increasing the DD velocity 
dispersion and decreasing the density we consider two further benchmark models, one with 
odd = 50kms _1 and pdd = 0.15/Oh and one with ctdd = 100 kms -1 and pr>D = 0.15/3H- In 
both cases we keep v\ ag = 50kms _1 and maintain the assumption of an isotropic gaussian 
velocity distribution. Ref. [4C] argues that the DD velocity distribution is better fit by a 
Tsallis distribution [74]. Given the substantial uncertainties in the DD density and velocity 
dispersion we do not investigate the effect of the uncertainty in the shape of the dark disc 
velocity distribution. For all three benchmark dark disc models, for simplicity and following 
Ref. 10711 ) we use the standard halo model with standard parameters for the dark matter 
halo and fix the total local density to the standard value pu + pdd = 0.3GeVcm -3 . 



3.3 Summary 

The benchmark models we consider are summarised in Table 1. For compactness we refer 
to the standard halo model with v c = 220 km s^ 1 and v esc = 544kms _1 as the standard 
halo model with standard parameters (SHSP). The normalised speed distribution for each 
of the models is plotted in fig. |l|. For the standard halo model, increasing (decreasing) v c 
increases (decreases) both vo, the value of v at which f{v) peaks, and also the width of f{v). 
The speed distribution of a pure DD has the same qualitative shape as the standard halo 
model. The effect of a DD on the total normalised speed distribution depends strongly on 
both the DD density and speed dispersion. With a high DD density and a low DD velocity 
dispersion there is a large additional peak in the speed distribution as low speed. As the 
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Figure 1: The normalised speed distributions. Top panel, the standard halo model: SHSP (solid 
line), SHu 0SC H (dotted), SHi; c L (short dashed) and SHt> c H (long dashed). Bottom left, modi- 
fied Maxwellian fits to VL2 simulation: SHSP (solid line) SIMsh (long dashed), SIMspmed (short 
dashed), SIMspl6 and SIMsp84 (dotted). Bottom right, dark disc models: SHSP (solid line), 
DDpHcrL (dotted), DDpLerL (short dashed), DDpLcrH (long dashed). Note the different scale in 
this and subsequent figures for the dark disc models. 



DD density is decreased the height of the peak decreases. If the DD velocity dispersion 
is increased the separation of the DD speed peak and halo speed peak decreases. For 
DDpLcrH, which has a low DD density and a large speed dispersion, the speed distribution 



has a single peak, at a lower speed than the standard halo model. As previously found [34] 
the speed distributions from the modified Maxwellian fits to the VL2 simulation data have 
less low speed and more high speed particles than the standard halo model with the same 
peak speed, vq. However the differences, in the lab frame, are fairly small p4] ]. The best 
fit to the shell data and the median fit to the sphere data are fairly similar. The scatter 
between sphere fits is 0(10%). If the simulation fits were compared to a standard halo 
with the same circular speed, v c (which is the observable quantity) rather than the same 
peak speed, vq, the deviations from the standard halo model would, however, be larger. 

4. Results 

4.1 Differential event rate 

In fig. ^ we plot T(f m ; n ), the time averaged value of the model independent parameterisation 



- 8 - 




of the differential event rate, eq.( |2.3| ), for each of the benchmark velocity distributions. For 
the standard halo model T(v m - m ) is approximately given by |nj 



T(v \ ~ n 220 kmS 1 PYn ( n V ^ \ (Mi 
1 [V m - m ) Ri a\ exp I — CL2 — 5- , (4.1 J 



where oi and a2 are constants of order unity i.e. increasing v c decreases both the overall 
normalisation and the rate at which T(v m ; n ) decreases with increasing v m { n . Varying the 
escape speed has a negligible effect, apart from at large v m i n . With a DD T(0) is larger, 
and the initial decrease in T(v m in) as v m i n is increased is more rapid. The smaller the DD 
density and the closer the speed dispersion to that of the standard halo, the smaller the 
difference from the standard halo. For the two models with a low DD density, the change 
in T(v m i n ) is relatively small for v m \ n > 100 km s -1 . Consequently in these cases there 
will only be a significant change in the mean differential event rate if the WIMP mass 
and/or experimental energy threshold are sufficiently low. For the modified Maxwellian 
fits to the VL2 simulation data the shape of T(v m { n ) is qualitatively similar to that of the 
standard halo model, but T(0) is smaller and the fall off of T(v m m) with increasing v m \ n is 
slower. The shell/median sphere fit are similar and the difference between them and the 
standard halo is similar to the spread in the sphere fits. The differences are small, however, 
compared with those from the uncertainty in the value of v c . We have checked that using 
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m ( (GeV) m x (GeV) 

Figure 3: Exclusion limits for an ideal experiment with a Ge target and an exposure £ = 10 3 kg day 
which detects zero events. Line types as fig. 1. 



the tabulated data from the VL2 simulation from Ref. [62] produces a similar uncertainty 
in T(u m i n ) as the spread in the sphere fits. 

In the absence of a positive signal experiments place exclusion limits on the WIMP 
mass and cross-section, which depend on the WIMP distribution |75|, pH] . To illustrate the 
effect of the uncertainty in the WIMP velocity distribution on exclusion limits we calculate 
exclusion limits for an ideal experiment (perfect energy resolution, perfect efficiency, zero 
background and zero energy threshold) using a Ge target with an exposure £ = 10 3 kg day. 
For each mass we find the cross-section for which the expected number of events is equal to 
three (the 95% upper confidence limit when zero events are observed) |76]| . The resulting 
exclusion limits are shown in fig. |3| for each velocity distribution. The total WIMP flux is 
proportional to the mean WIMP speed, which for the standard halo model is proportional 
to v c . Therefore the naive expectation is that increasing v c increases the total event rate 
and hence leads to a tighter constraint on a p . As can be seen in fig. || this is the case for 
small m x . For larger m x the, energy dependent, suppression of the event rate by the form 
factor, means that when v c is increased the total event rate in fact decreases and the limit 
on dp becomes weaker (see also Ref. |5^1 ). Increasing the energy threshold weakens the 
constraints, and increases the mass at which the transition occurs. Changing the escape 
speed only affects the exclusion limits significantly for light (m x < O(10 GeV)) WIMPs, 
see Ref. [41|. For a high density DD the increase in T(v m i n ) for small u m ; n means that the 
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Figure 4: Fractional mass limits for an ideal Ge detector and an exposure 8 = 3 x 10 4 kgday 
calculated assuming the standard halo model with standard parameters. Line types as fig. 1. 



exclusion limit is a factor of a few weaker for all m x . For the low density DD models and 
the simulation fits the change in the exclusion limit is, as expected, relatively small. 

Energy dependent efficiency and/or background events will change the shape of the 
exclusion limits, see e.g. Ref. |75|, fjl[ . As discussed above the uncertainty in the local 
density translates directly into an uncertainty in a p which is the same for all experiments. 
Note that, as pointed out in Ref. |4i| , the values of p x and v c are correlated. 

Once events are detected the shape of the energy spectrum can be used to measure 
the WIMP mass [|ll], 12, [l3|, [b|, ff5|, [R|. Uncertainty in the velocity distribution leads to a 



systematic uncertainty in the determination of the WIMP mass [14]. We use the method 



used in Refs. [14, 16] to calculate the limits on the WIMP mass which would be obtained 
for each velocity distribution, if the data were analysed assuming the standard halo model 
with standard parameters. As before we assume an ideal Ge detector 2 a WIMP cross- 
section, Up = 10~ 8 pb and an exposure of £ = 3 x 10 4 kg day. We estimate the WIMP mass 
and cross-section by maximising the extended likelihood function, e.g. Ref. [ffTj: 

t= *y »irw). (4.2) 



2 See Ref. ]lq | for an exploration of how varying the detector capabilities affects the WIMP mass deter- 
mination. 



Here iV ex p t is the number of events observed, E\{i = 1, N expi ) are the energies of the 
events observed, f(E) is the normalised differential event rate and A = £ J^°(dR/dE)dE 
is the mean number of events. We calculate the probability distribution of the maximum 
likelihood estimator of the WIMP mass, for each input WIMP mass, by simulating 10 ex- 
periments and finding the 2.5%, 50% and 97.5% percentiles, m^ m , of the mass distribution. 

The fractional mass limits, (m^ m — m™)/m™ are plotted in fig. || as a function of the 
input WIMP mass, m™. As discussed in Refs. 0, [14], @, assuming an erroneous value of 



Am Y 

± 



1 + 



(4.3) 



v c leads to a systematic error in the mass determination which increases with increasing 
m x : 

m A J\ v c 

A dark disc has a similar effect. There is a population of WIMPs with lower speeds than 
assumed, and hence the WIMP mass is systematically underestimated. The systematic 
underestimate is substantial (~ 10 — 50%, increasing with increasing m x ) if the DD density 
is large and the speed dispersion is significantly different from that of the halo. For the 
two models with a low DD density, the systematic underestimate is smaller, ~ 5 — 20%. 
The larger width of the modified Maxwellian speed distribution leads to a systematic 
overestimate of the WIMP mass in the range ~ 2 — 10%, increasing weakly with increasing 
WIMP mass. 



4.2 Annual modulation 

We consider the amplitude of the modulation in the model independent parameterisation 
of the differential event rate, AT(v m - m ) = m&x{T(v m - m ,t) — T(t> m i n )], 3 and the date t p on 
which this maximum occurs. For small f m i n , the maximum event rate occurs in Winter ff"8|j . 
As 

^min is increased, AT(t; m i n ) initially decreases to zero at which point the phase of 
the annual modulation changes rapidly and the maximum occurs in Summer. As v m i n is 
increased further the amplitude increases one more to a local maximum (which we refer to 
as the Summer maximum), before decreasing again and tending to zero ffq , [79| ]. 

In figs. [| and ^ we plot AT(v m ; n ) and the day of the year on which the maximum occurs, 
t p , respectively, for each of the benchmark velocity distributions. We do not include a plot 
of t p for the standard halo model as the low and high v m { n values of t p only change by ~ 1 
day as v c is varied. 

For the standard halo model, as v c is increased the most likely velocity and the width of 
the velocity distribution both increase. Consequently AT(0) and the Summer maximum of 
AT(v m i n ) both decrease, the value of v m i n at which the Summer maximum occurs increases 
and the decline in the amplitude for large v m i n is less rapid. The value of v m { n at which 
the maximum switches from Winter to Summer increases. 

For a pure DD (no halo) the behaviour would be qualitatively similar to the standard 
halo model. Due to the speed distribution peaking at a smaller speed, and having smaller 
dispersion, AT(0) and the Summer maximum would be larger, and the phase change would 
happen at smaller f m i n , ~ 60kms _1 . When a DD is added to the standard halo the net 

3 Some papers, e.g. Ref. |34], consider the fractional annual modulation, which is largest for large u m m. 
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effect is more complicated. For DDpHcL there are two local Summer maxima in the 
amplitude, a large one at v m - m ~ 110 km s -1 from the DD and a smaller one from the halo 
in the usual position, w m ; n ~ 330 km s -1 . In between 80kms _1 < w m ; n < 200 km s" 1 the 
maximum occurs at t p ~ 125 days. For DDpLcrL the DD density is not high enough to 
produce a second Summer maximum, however is does lead to a significant variation of t p for 
80kms _i < 

^min ^5 200 km s . This is because for these values of v m i n the contributions 
from the DD and the halo are out of phase, and partly cancel. For DDpLuH the change 
from the standard halo model is fairly small. 

For the modified Maxwellian distributions the qualitative behaviour of AT(v m ; n ) is 
the same as for the standard halo model; AT(0) is smaller, the Summer maximum occurs 
at a smaller v m i n , but its amplitude can be either smaller or larger. The changes in the 
amplitude of the modulation are larger than those in the mean differential event rate. 
However, the changes from the uncertainty in v c are again larger than those from the 
uncertainty in the shape of the speed distribution (unless there is a dark disc with small 
speed dispersion). The change in t p from Winter to Summer happens more gradually, and 
for large v m ; n t p is typically a few days smaller. For both the SHSP and the modified 
Maxwellian fits AT(~ 200 km s -1 ) ~ and hence the corresponding value of t p is hard to 
calculate accurately, and not particularly physically meaningful. 
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Figure 6: Day of the year, t p , on which the maximum of the differential event rate, AT(u m ; n ) = 
max[T(w min , t) — T(w min )], occurs. Line types as Fig. f for the modified Maxwellian fits to the 
simulation data (left panel) and dark disc models (right panel). 
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Figure 7: Contour plots of the model independent parameterisation of the direction dependence, 
T(fmin, cos 7), for the standard halo model. Top row: SHSP and SHu 0SC H, bottom row: SHu c L and 
SHi; c H. Contours have spacing 0.5 between and 3.5. 



4.3 Direction dependence 

The model independent parameterisation of the direction dependence, T(f m in, cos 7) as 
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defined in eq. ( |2.9| ), is shown in figs. |?], || and ^ for the standard halo model, the modified 
Maxwellian fits to the simulation data and the dark disc models respectively. For the 
standard halo model 101 



T(f m in,cos7) « 2y/nexp 



orb,p 



+ Vq) cos 7 ■ 



(4.4) 



where t>e rb,P is the component of the Earth's velocity parallel to the direction of Solar 
motion. This has a maximum value, for v m \ n = (i>° rb ' p + W0)cos7, of T ~ ~ 3.5. 

As can be seen in fig. |7], for the standard halo model T(u m i n , cos 7) is constant for fixed 
Wmin/ cos 7. This is not the case for the other velocity distributions. For the anisotropic 
modified Maxwellian distribution, around the peak of the T{v m \ n , cos 7) distribution as v m \ n 
is decreased, with v m \ n / cos 7 fixed, the value of T decreases, while the contours of low fixed 
T are convex. For the DD models the additional population of low speed WIMPs means 
that T peaks at smaller v m { n and has a larger peak value (~ 6.5 for DDpHerL compared 
with w 3.5 for the standard halo). Around the peak of the T{v m \ n , cos 7) distribution in 

this CclSG clS ^miri IS 

decreased, with v m \ n / cos 7 fixed, the value of T decreases. As before, 
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Figure 9: As fig. [7] for the dark disc models (and, for comparison, the standard halo model with 
standard parameters). Top row: SHSP and DDpHerL, bottom row: DDpLuL and DDpLerH. For 
clarity these plots have a different scale to figs. ^ and ^. Contours have spacing 0.5 between and 
6.5. 



the smaller the DD density and the closer the speed dispersion to that of the standard 
halo, the smaller the difference from the standard halo. 

We use the methods presented in Ref. [29| to calculate the number of events required 
to determine that the recoil distribution is not isotropic if m x = 100 GeV. Briefly, we use 
the HADES code [ 30 1 to simulate the recoil direction distribution for a S detector with 
energy threshold 20keV 4 . We assume that the recoil directions, including their senses, are 
reconstructed perfectly in 3d and the background is zero. These are optimistic assumptions 
and therefore our results provide a lower limit on the number of events required by a real 
detector. For 3-d data the most powerful test for rejecting isotropy uses the average of the 
cosine of the angle between the direction of solar motion and the recoil direction, (cos 7) [29]. 
We calculate the probability distribution of (cos 7), for a given number of events N, by 
Monte Carlo generating 10 4 experiments and find the number of events required to reject 
isotropy at 95% confidence in 95% of experiments, A r i so . For further details see Ref. p9fl. 



4 It is impossible for a directional detector to have zero energy threshold, even in principle, as the lengths 
of the nuclear recoils tend to zero in this limit, and hence their direction is unmeasurable. 



For all the velocity distributions considered iV iso = 9. 

The WIMP origin of an anisotropic recoil distribution could be checked by measuring 



the median recoil direction [81|. For a smooth WIMP distribution the median inverse recoil 
direction coincides with the direction of solar motion, modulo statistical fluctuations. The 
median inverse direction is defined as the direction x me d which minimises the sum of the 
arc- lengths between x me d and the individual inverse recoil directions x; [82]. It is found by 
minimising 



N 



M = 5Z C ° S l ( X med-Xj 



(4.5) 



where N is the number of events. We determine the number of events required to confirm 
the direction of solar motion as the median inverse recoil direction at 95% confidence using 
the distribution of A, the angle between the median direction and the direction of solar 
motion, x Q : 

A = cos~ 1 (x mcd .x Q ) . (4.6) 

We calculate the probability distribution of A, for a given number of events N, by Monte 
Carlo generating 10 4 experiments and find the number of events required to reject the 
median direction being random at 95% confidence in 95% of experiments, N me ^. For 



further details see Ref. | 8i| . The value of N me< ±, is only weakly dependent on the velocity 
distribution. For SHv c H, N me d = 32 while for all the other velocity distributions N me( i = 27 
or 28. 

The limited variation of N[ so and N me ^ show that the directional signals are robust to 
(plausible) uncertainties in the velocity distribution. The modified Maxwellian fits to the 
simulation data do not include the stochastic features found at high speeds however. These 
features can cause the median inverse direction of high energy recoils to deviate from the 



direction of solar motion [34]. As discussed in Ref. [Bl| the deviation will be small unless 
the detector is only sensitive to the high speed tail of the speed distribution (i.e. if the 
WIMP mass is small and/or the energy threshold is high). Studying the median direction 
of high energy recoils could in fact allow high speed features to be detected, hence probing 
the formation history of the Milky Way. 



5. Summary 

Direct detection event rate calculations often assume the standard halo model, with an 
isotropic Maxwellian velocity distribution. It is well known, however, that the energy [17, 
1|, [75|, ||, |37|, |2|, H§ 



1| , time H, |g, H |2|, |§ ||, g|, g|, H> HI' HI H 



40] 



and direction [27, Bq, 2£, 34] dependence of the direct detection event rate all depend on 



the local WIMP distribution. We have updated these studies in light of recent numerical 
simulations and observational measurements of the local circular speed and escape speed. 

The local circular speed has a model dependent relation to the radial velocity dis- 
persion. Recent determinations, e.g. Refs 51, |53|| , have relatively small statistical errors, 
however Ref. 54] found values ranging from v c = 200 to 280kms _1 , depending on the 
model of the MW assumed. The most recent determination of the local escape speed [55] 
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finds v esc = 544 km s" 1 , slightly lower than the historical standard value, v esc = 650 km s -1 . 
Similarly for a given model for the MW it is possible to determine the local DM density, 
p x , to ~ 10% accuracy |5(], j57|1 , however the systematic errors are likely to be substantially 
larger §|, [59], 

High resolution dark matter only simulations of the formation of Milky Way like dark 
matter halos, typically find speed distributions which can be fit with a modified Maxwellian 
distribution which is broader than the standard Maxwellian distribution |3(], [32], |34fl . 
Some simulations which include baryonic physics find that late merging sub-halos are pref- 
erentially dragged towards the disc, where they are destroyed leading to the formation of a 



co-rotating dark disc [31, |33|, 35]. The significance of a DD for direct detection experiments 



depends on its density and velocity distribution, which are highly uncertain [[37], 73, 40 1 . 

We have considered three types of velocity distribution: standard Maxwellian, mod- 
ified Maxwellian, standard Maxwellian plus dark disc. In each case we use a range of 
parameter values motivated by recent observations and simulations (as discussed in detail 
in Sec.|3l). 



Differential event rate and exclusion limits 

The systematic uncertainty in the local circular speed, v c , leads to a 0(10%) uncer- 
tainty in the differential event rate, and hence exclusion limits. A high density DD also 
produces significant changes in these quantities. A low density DD will only have a signifi- 
cant effect if the WIMP mass and/or energy threshold are sufficiently low. The dependence 
on the detailed shape of the velocity distribution is small. 

The normalisation of the differential event rate is directly proportional to the product 
of the local density and the WIMP cross-section. Therefore the uncertainty in the mea- 
surement of the local density propagates directly into an uncertainty on measurements of, 
or constrains on, the cross-section. 



Mass determination 

We studied the systematic errors in determinations of the WIMP mass which would 
occur if data from a future SuperCDMS like experiment is, erroneously, analysed assum- 
ing the standard halo model with standard parameters. Assuming an incorrect value of 
v c leads to a systematic error in the mass determination which increases with increasing 



in 



x. 



[13, 14, [16|. With a DD there is a population of WIMPs with lower speeds than as- 



sumed and hence the WIMP mass is underestimated. The size of the systematic error 
varies, with increasing m x , from 10 — 50% for a high density dark disc and from 2 — 10% 
for a low density dark disc. The larger width of the modified Maxwellian distribution leads 
to a smaller overestimate of m v . 



Annual modulation 

The annual modulation is more sensitive to the velocity distribution than the mean 
differential rate. Changing the value of v c or a DD with small speed dispersion can change 
the amplitude by a factor of order unity, while the uncertainty in the shape of the halo 
velocity distribution changes the amplitude by ~ O(10%). 
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The phase of the modulation depends only weakly on v c and the shape of the halo 
velocity distribution. With a small speed dispersion DD the phase at moderate energies 
changes by ~ 10 (100) days if the DD density is low (high). 

Direction dependence 

The detailed direction dependence of the event rate is sensitive to the velocity distri- 
bution, however the directional signals are robust. The number of events required to detect 
anisotropy does not change, while the number of events required to demonstrate that the 
median inverse recoil direction coincides with the direction of solar motion varies by of 
order 10%. 



Even with recent improvements in the resolution of numerical simulations and observa- 
tional determinations of the dark matter parameters, there are still significant uncertainties 
in the direct detection signals and WIMP parameter determinations. In particular the exis- 
tence of a DD could have a significant effect, depending on its density and speed dispersion. 

Considering a range of, data motivated, benchmark models is an improvement on sim- 
ply assuming the 'standard halo model'. However it is still not completely satisfactory. In 
particular it does not provide a formal analysis of errors. Several approaches to dealing 
with the impact of astrophysical uncertainties on direct detection data analysis have re- 



cently been proposed. Strigari and Trotta [83] have suggested using astronomical data and 
a model for the Milky Way mass distribution in a Monte Carlo Markov Chain analysis of 
direct detection data. Peter ]84j] has presented an approach which involves combining data 
sets from different direct detection experiments and jointly constraining a parametrisa- 
tion of the WIMP speed distribution and the WIMP parameters (mass and cross-section). 
These are promising directions, however in both cases the current implementations assume 
an overly restrictive form for the velocity distribution (a single isotropic Maxwellian). 
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